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Single-molecule magnets weakly coupled to two ferromagnetic leads act as memory devices in 
electronic circuits — their response depends on history, not just on the instantaneous applied voltage. 
We show that magnetic anisotropy introduces a wide separation of time scales between fast and 
slow relaxation processes in the system, which leads to a pronounced memory dependence in a wide 
intermediate time regime. We study the response to a harmonically varying bias voltage from slow 
to rapid driving within a master-equation approach. The system is not purely memristive but shows 
a partially capacitive response on short time scales. In the intermediate time regime, the molecular 
spin can be used as the state variable in a two-terminal molecular memory device. 

PACS numbers: 73.63.-b, 85.35.-p, 73.23.Hk, 05.60.Gg 



I. INTRODUCTION 

Electronic transport through magnetic molecules 
has been extensively studied both experimentally^^ 
and theoretically!^"^' xhe magnetic moment of these 
molecules can be realized by one or more transition-metal 
ions with organic ligands, by organic radicals, or by en- 
dohedral fullerenes. Transition-metal ions in a molecular 
ligand field typically show sizable magnetic anisotropy. 
In the case of a pronounced easy axis, one then speaks of 
single-molecule magnets (SMMs). For memory devices, 
easy-axis anisotropy is desirable since it introduces an 
energy barrier for spin reversal and thereby stabilizes the 
spin in the up or down orientation.^ 

Driving the system with an external electric field al- 
lows the control and manipulation of the molecular spin 
state, and consequently writing and reading of informa- 
tion using SMMs.^ Clearly, the response of the system 
to the applied electric field is not purely instantaneous — 
the system has memory. This can be compared to the 
response of SMMs to a magnetic field. At low tempera- 
tures, they show pronounced hysteresis in the magneti- 
zation (of bulk crystals of noninteracting SMMs) vs the 
magnetic field. 1 35 1 36 * 

Similarly, if a harmonically varying bias voltage is ap- 
plied, we expect hysteresis loops in the spin, charge, 
and current dynamics vs the applied voltage to de- 
velop, whose amplitude depends on the voltage ampli- 
tude and frequency — as in other systems whose spin 
polarization can be controlled by a voltage or current; 
see, e.g., Ref. 37. At very low frequencies, the spin 
dynamics can easily follow the external field and lit- 
tle or no hysteresis is expected. At very high frequen- 
cies, the spin dynamics are "frozen." There is, how- 
ever, an intermediate frequency range — comparable to 
the inverse spin-relaxation time(s) — where the hysteresis 
is most pronounced. In this range, the device state, at 
any given time, is strongly dependent on the history of 
states through which the system has evolved. In that 



case, we expect the resistance of the device to be a func- 
tion of the state variable x that describes its memory 
(the spin polarization) and possibly of the protocol with 
which the voltage V(t) has been applied, i.e., the en- 
tire wave form of V(t). In other words, the resistive 
response can be characterized by a function of the type 
R(x, V, t). Such a device goes under the name of memris- 
tive (for "memory resistive") systemPsJSQl Resistors with 
memory are experiencing a surge of research activity, in 
part due to promising applications in memory storage, 
but also because of their ubiquity in diverse areas rang- 
ing from nontraditional computing to biophysics It 
is natural to think that SMMs form another example 
of memristive systems, with the molecular spin play- 
ing the role of internal state variable. This would have 
the added advantag e of c ombining memristive and spin- 
tronics functionality^^ in molecular junctions. Indeed, 
Miyamachi et al. 13 have recently demonstrated memris- 
tive behavior of single Fe(l, 10 — phcnanthrolinc)2(NCS)2 
molecules on CuN under a scanning tunneling micro- 
scope. In their case, the molecule is switched between 
a high-spin (S — 2) and a low-spin (S — 0) state of the 
Fe 2+ ion, which is connected with a change in conforma- 
tion and conductance. Our case is quite different: We 
consider the switching of a molecular local spin of fixed 
length S over an anisotropy barrier. 

In this paper, we show that the response of this sys- 
tem is only partially memristive. In addition, capacitive 
components emerge on short time scales. The capaci- 
tive components are related to the charging energy of the 
molecule and thus to the Coulomb repulsion of electrons. 

Transport through magnetic molecules is typically 
dominated by a strong exchange interaction between the 
spin of mobile electrons and the local molecular spin, in 
addition to the large Coulomb repulsion. We are thus 
faced with a strongly interacting nonequilibrium system, 
which makes a quantitative description difficult in gen- 
eral. However, if the coupling between the molecule and 
the metallic leads is weak, as is often the case in break 
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junctions, this coupling can be used as the small para- 
meter in a perturbative approach. This can be done in 
the framework of the master equation, which has the ad- 
vantage that the strong interactions within the molecule 
can be treated exactly. The master equation has been 
applied to t ransport throug h magnetic molecules by var- 
ious groups i^- 4 * 1 5 * 1 7 H HEimi] It provides an ensemble de- 
scription, not a description of individual time series of 
single-molecule devices. Statements we make about the 
memristive properties are thus to be understood in an 
ensemble sense. On the other hand, devices consisting 
of a monolayer of weakly interacting molecules between 
metallic electrodes, in which many molecules conduct in 
parallel, are self-averaging. In this case, we predict the 
time-dependent observables of a single device. 

The analysis of SMMs as memory devices requires us 
to study their dynamics under a time-dependent bias. 
In Refs. [17] and [HI their relaxation for constant or sud- 
denly switched bias has been considered. Here, we con- 
sider the current, charge, and spin response to a har- 
monically varying bias V(t) = Vq sinwi, which is easily 
realizable experimentally. For a non-magnetic molecule 
involving a vibrational mode, the response to a harmonic 
voltage has been recently studied by Donarini et al^ In 
this case, Ftanck- Condon blockade leads to interesting 
dynamics ! 45 l 46 l 

The remainder of this paper is organized as follows: In 
Sec. [Tl] we present our model, followed by a discussion 
of the master-equation approach in Sec. |III| The results 
are presented and discussed in Sec. [TV] Finally, in Sec.[V| 
we summarize the main points, address possible limiting 
effects, and draw some conclusions. 



II. MODEL 

Our device consists of a magnetic molecule coupled to 
two ferromagnetic leads. The full system is described by 
the Hamiltonian H = H mo \ + H \ ea g s + -ffhybj where the 
molecular Hamiltonian read^ilLL^MI 



e d ^ d ^d a 



Ud\d t dld i -Js-S-K 2 {S z )' z , (1) 



where d\ (d a ) creates (annihilates) an electron of 
spin er in the molecular orbital with energy e^, s = 
J2aa' dt-((Tacr> /2)d a is the corresponding spin operator, 
and S is the spin operator of a local spin of length 5. 
The two spins are coupled by the exchange interaction J 
and the local spin is subject to an easy-axis anisotropy 
of strength K% > 0. The extension of this model by in- 
cluding more than one electronic orbital, or more than 
one local spin, does not pose any conceptual difficulties. 
In a break-junction setup, the on-site energy e c i could be 
tuned by a gate voltage. However, this tunability is not 
necessary for our conclusions. 

The molecular Hamiltonian H mo \ commutes with the 
z component of the total spin S tot = s + S so that the 
eigenvalue m of 5 t z ot is a good quantum number. We 
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FIG. 1: (Color online) Molecular energy levels -E+m vs mag- 
netic quantum number m for ed = 0.2 eV, U = 10 eV, J = 
0.1 eV, Ki = 40meV, and S = 2. Levels for electron occupa- 
tion number n = (n = 1) are shown as heavy blue (medium 
red) bars. Levels for electron number n — 2 are outside 
the range of the plot. The arrows indicate the lowest-energy 
sequential-tunneling transitions out of the ground states. 



show in Fig. [TJ the energy levels of -ff m oi vs m for the 
parameter values td = 0.2 eV, U = 10 eV, J = 0.1 eV, 
K2 — 40meV, and S = 2, which we will also use below 
to illustrate our results. Note that we are using a large 
value of the anisotropy energy in order to most clearly 
show the generic behavior. For the smaller anisotropies 
of SMMsp^the interesting physics would occur in a nar- 
rower voltage range and at lower temperatures. For n = 
and for n — 2 electrons on the molecule, the total spin is 
just given by the local spin S so that there are 2S+1 levels 
in these charge sectors. For one electron, n = 1, its spin 
1/2 combines with the local spin S to form two multiplets 
of 2(5 -1/2) + 1 = 25 and 2(5 + 1/2) + 1 = 25 + 2 states. 
The splitting between the two multiplets is on the order 
of JS. The easy-axis anisotropy leads to the parabolic 
dispersion seen for all multiplets in Fig. [TJ The total di- 
mension of the molecular Fock space is Np — 4(25 + 1). 
The leads are described by the Hamiltonian 



leads — / 



(2) 



where c^, ko . (c„k<r) creates (annihilates) an electron with 
wave vector k and spin a in lead v — L, R. Below, we will 
assume that the leads are ferromagnetic with opposite 
magnetizations parallel to the z direction. The molecule 
and the leads are coupled by the hybridization term 



hyb 



^hyb 



(3) 



where, for simplicity, we have assumed the tunneling ma- 
trix element thyb to be real and independent of the lead 
index, wave vector, and spin. The number of sites, N, 
in each lead is introduced by the Fourier transformation 
into momentum space and drops out of the physical re- 
sults. 
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III. MASTER EQUATION 

The master-equation approach starts from the exact 
von Neumann equation (we set H = 1), 



dp 
~dt 



= -i[H,p] 



(4) 



for the density (statistical) operator p for the complete 
system. The reduced density operator of the molecule 
is obtained by tracing out the degrees of freedom of the 
leads, jO mo i = Triads P- The resulting equation of motion 
is the master equation for p mo \W^^ln principle, an exact 
master equation that is local in time, 



dt 



Pmol 



(t)]-K(t,to) Pmol (t), (5) 



can be derived even for time-dependent Hamiltonians if 
the molecule and the leads were in a product state at 
some initial time to The first term on the right-hand 
side of Eq. ^ describes the time evolution of the decou- 
pled molecule, whereas the second term involving a linear 
superoperator 1Z describes relaxation. 

In practice, approximations are needed to obtain su- 
peroperators 1Z that preserve the positivity of the density 
operator. Here we make three simplifications: (i) We em- 
ploy the sequential-tunneling approximation, which con- 
sists of keeping only terms up to second order in the tun- 
neling matrix element ihyb- This approximation is rea- 
sonable if T <C fcflS, where O is the temperature, (ii) We 
only consider reduced density operators p mo \ that are di- 
agonal in the eigenbasis of H mo \. Actually, this is not an 
approximation: Below, we will always assume that the 
time evolution starts from equilibrium or from a pure 
diagonal state. If we start from such a diagonal p mo \, 
the exact time evolution does not generate nonzero off- 
diagonal components. This is because off-diagonal com- 
ponents would correspond to superpositions (coherences) 
of states with different charge or different spin component 
St ot . The former would break U(l) charge symmetry, 
and could only be expected in superconducting systems. 
The latter would lead to nonzero averages of S* ot or S% ot , 
which would require spontaneous breaking of the spin- 
rotation symmetry around the z axis, (iii) We employ 
the Markov approximation, which posits that correlation 
functions describing the memory of the leads decay on a 
time scale ri ea ds that is much shorter than all experimen- 
tally relevant time scales. This requires the oscillation 
period T of the bias voltage to satisfy T ^> ri ca d s . Since 
the leads are metals with typical relaxation times in the 
femtosecond range, this is easily fulfilled. The Markov 
approximation is necessary here since it allows us to use 
the instantaneous value of the bias voltage in the master 
equation. 

The derivation of the resulting master equation is 
standar d 14 * 17 * 46 * 54 ^^ and we only present the results. We 
assume the bias voltage V to be split evenly between 
the two molecule-lead contacts. The two leads are as- 
sumed to be identical ferromagnetic metals with their 



magnetizations along the z direction but of opposite 
sign. The relevant parameter is the ratio p = -/V m ; n /./V lna j 
between the densities of states (assumed to be energy 
independent) for minority-spin and majority-spin elec- 
trons. We write the diagonal reduced density operator 
as p mo i = diag(Pi, P%, ■ ■ • ), where the Pi are the prob- 
abilities of many-particle eigenstates |i) of H mo \. The 
master equation then takes the form 



dP 
dt 



(6) 



with Rj->i the transition rate from state \j) to state \i) 
due to sequential tunneling. The rate can be written as a 
sum over contributions from spin up and spin down and 
from the left and right leads, 
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(7) 
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(10) 



(11) 



where f(E) is the Fcrmi-Dirac distribution function, Ei 
is the eigenenergy of molecular state \i), Dfj = (i\d a \j) 
are matrix elements of the electron annihilation operator 
between molecular eigenstates, and T = 27r|th y b| 2 -^maj 
quantifies the coupling of majority electrons to the leads. 

The average occupation number, the z component of 
the electron spin, and the z component of the local spin 
are given by 



^Piiil^dtd^i) 



(s z ) = E p * 



dldf — d\d\ 



(12) 

(13) 
(14) 
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respectively. The average charge current between lead v 
and the molecule 



(15) 



where the numerical value of v is +1 (—1) for the left 
(right) lead and rii = (i\ J2 a °%d a \i). The current is 
counted as positive if it is flowing from left to right. 

While for the stationary state the left and right cur- 
rents are equal, this is not generally the case for time- 
dependent Pinoi-^ It is then crucial to realize that an 
ammeter in, say, the left lead nevertheless measures the 
symmetrized current, 



/ : 



(II) + (In) 



(16) 



The reason is the following: (I v ) only contains the tun- 
neling (particle) current through the contact between the 
molecule and lead v. In addition, there are displacement 
currents across the contacts resu lting from charging of 
the molecule-lead capacitors! 57 * 58 ! The displacement cur- 
rents are equal to the charging currents, as seen from the 
simple example of a pure capacitor. An ammeter placed 
in the left lead picks up the sum of the tunneling and the 
charging currents. By recalling that for a symmetric de- 
vice the displacement currents in both barriers are equal 
in magnitude but opposite in signp^J and that the sum of 
particle and displacement current is divergence free, one 
can show that the sum of the particle and displacement 
currents in the left contact, and thus the current in the 
ammeter, equals the symmetrized particle current!^ 

The master equation ^ is solved numerically by dis- 
cretizing time and propagating the probabilities Pi(t) for- 
ward step by step. We will compare our results to the 
stationary solution at constant bias voltage, which consti- 
tutes the limit of infinitely slow driving. The stationary 
solution P°° is found by setting the left-hand side of Eq. 
|6]) to zero, resulting in 







(17) 



with 



A, 



-Ei^-Ri-Hfe for i= i- 



(18) 



The stationary probability vector P°° = (Pj 30 , P£° , . . .) 
is thus the right eigenvector of the transition-rate matrix 
A with zero eigenvalue. A has at least one vanishing 
eigenvalue, i.e., one stationary state, since (1,1,..., 1) is 
clearly a left eigenvector with zero eigenvalue. 

We can also conclude that since our model is ergodic — 
any state \i) can be reached from any other state by a fi- 
nite number of sequential-tunneling transitions with non- 
vanishing rates — the stationary state is unique. However, 
the matrix A is often ill conditioned since the rates Rj-ii 



are spread over many orders of magnitude due to the 
Fermi functions. This leads to the problem that diago- 
nalization routines working with machine precision ob- 
tain more than one eigenvalue that is numerically zero, 
and thereby fail to find the unique stationary state. We 
overcome this difficulty by diagonalizing A with high pre- 
cision in mathematical We use LU decomposition to 
estimate the L°° condition number and adapt the number 
of digits kept in the diagonalization so that the resulting 
P°° contains at least 12 significant digits. 

For periodic bias voltage V(t), the system will not re- 
lax toward a stationary state but will approach a peri- 
odic cycle. We define the time-evolution matrix II for 
one full period by P(t + T) = II P(t). To make II unique, 
we choose the time t in such a way that V(t) = and 
V'(t) > 0, i.e., at vanishing bias voltage on the upsweep. 
Il is the stochastic matrix of a discrete-time Markov pro- 
cess. The periodic state is characterized by the prob- 
ability vector pp cr = (P 1 pcr , P|' cr , . . .) which is a right 
eigenvector of the stochastic matrix with eigenvalue one, 
ppcr _ jrpP 01 '. The stochastic matrix II has at least one 
unity eigenvalue, and this eigenvalue is nondegenerate, 
since the discrete-time process is ergodic, in analogy to 
the discussion for A above. Starting from pp cr , the pe- 
riodic time-dependence can be obtained by integrating 
the master equation ([6]) over one period. The stochastic 
matrix II is evaluated numerically by discretizing time 
as II = no<t<T[-"- + AtA(i)], where the transition-rate 
matrix A now depends on time through the bias volt- 
age. We normalize II after each matrix multiplication 
such that ^2 m II m „ = 1 for all n, which is required to 
conserve probability. II can also be ill conditioned and 
we apply the method sketched above to find the unique 
periodic state. 



IV. RESULTS 

Before analyzing the dependence on the period and the 
amplitude of the bias voltage in detail, let us show the 
typical behavior of our system for one parameter set. Fig- 
ure[2]shows the approach to the periodic regime when the 
system is initialized in its equilibrium state at time t = 
and then driven by a bias V(t) — Vq s'mujt with period 
T = 2tt/uj. The current, occupation number, and the z 
component of the total spin approach periodic behavior 
within a few periods. As expected, all three observables 
show hysteresis. Since the state of the system evidently 
depends on its history, we immediately see that our sys- 
tem is indeed a memory device.^ Its internal state is 
described by the probability vector P(t), which contains 
Np — 1 independent real variables, where Np is the di- 
mension of the molecular Fock space. However, it is not 
a purely memristive system: A memristive system would 
satisfy equations of the forrrPSl 



J(t) = G(V(t),P(t))V(t), 



(19) 
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FIG. 2: (Color online) (a) Current, (b) occupation number, 
and (c) z component of the total spin vs bias voltage for 
the SMM driven by a harmonic bias voltage of amplitude 
Vo = 0.5 V and period T = 10 T" 1 . Ten full periods are 
shown as solid blue (gray) curves. The arrows indicate the 
directions in which the loops are traversed. The black dashed 
curves refer to the stationary state, corresponding to T — > oo. 
The spin polarization in the leads is p — N m i n /N ma j = 0.5. 
The other model parameters are ed = 0.2 eV, U = 10 eV, 
J = 0.1 eV, K 2 = 40meV, and S = 2, as in Fig. [T] 



dP 

lit 



F(V(t),P(t)). 



(20) 



In our case, the second equation is the master equation 
On the other hand, our Eqs. (15 1 and ( fl6| for the 
current cannot be written in the form of Eq. ( |19[ ) for all 
V. Equation ( 19 ) implies that the current vanishes for 
zero bias for a mcmristive system, whereas Fig. [2] shows 
that it does not vanish for our device — the hysteresis loop 
is not pinched at V = 0. Physically, this is because we 
have additional capacitive or inductive effects. The I- 
V characteristics of a harmonically driven pure capac- 
itor (inductor) show an ellipse that is traversed in the 



clockwise (counterclockwise) direction. Since the loop in 
Fig. [2^ a) is clockwise, the behavior is partially capaci- 
tive. This is reasonable since the molecule is transiently 
charged, as shown in Fig. [2|b). 

For comparison, Fig. [2] also shows the voltage depen- 
dences of all observables in the stationary state (dashed 
curves). Note that the average spin in the stationary 
state is nonzero and depends on the bias voltage solely 
because the leads are magnetically polarized. If electrons 
are predominantly moving from left to right (I < 0), the 
left lead injects predominantly spin-up electrons, while 
the right lead absorbs predominantly spin-down elec- 
trons. The result is a positive spin polarization on the 
molecule, as seen in Fig. [2jc) . The stationary curves all 
show plateaus separated by thermally broadened steps. 
The dynamical current and charge mainly relax toward 
the stationary values at the plateaus but cannot follow 
the rapid changes at the steps. The visibility of the re- 
laxation suggests that the driving period is comparable 
to the relevant relaxation times. On the other hand, the 
spin hysteresis loop bears no resemblance to the station- 
ary curve, showing that the spin cannot relax rapidly 
enough to approach its stationary value. We will discuss 
these points in what follows. 



A. Dependence on the voltage period 

To elucidate the role of the driving period T, we have 
determined the periodic behavior for various T's as de- 
scribed in Sec. |III[ keeping all other parameters fixed. 
Figures |3j |4j and [5] show single hysteresis loops for cur- 
rent, charge, and spin in the periodic regime. Since the 
hysteresis in the current is not very pronounced on the 
scale of the figures, we show an enlargement close to the 
voltage maximum in Fig. [6j The natural time scale is 
r _1 = (27r|ihyb| 2 ^maj) _1 , the inverse of the character- 
istic tunneling rate. Note that f or a typical current of 
200 pA in the transmitting regimep^the tunneling rate 
is on the order of V « 10 9 s _1 . In the limit of rapid 
driving (T — > 0), shown in Fig. |3j all hysteresis loops 
close since the molecule cannot follow the rapid bias. The 
charge and the spin approach constant values determined 
by the time-averaged rates, 



Rr, 



dt R n 



>(*)■ 



(21) 



The current loop also closes but it does not become a 
horizontal line since the current in Eq. (15 1 explicitly 



depends on the instantaneous rates, which in turn depend 
on the instantaneous bias. 

We note that the hysteresis loops for the current, in 
particular at rapid driving, show additional steplike fea- 
tures absent from the stationary curves. These features 
are due to excited-state-to-excited-state (ETE) transi- 
tions. Such transitions are also observed in the station- 
ary stated but only if the energy of the ETE transition 
is higher than the energy of the transition populating 
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FIG. 3: (Color online) (a) Current, (b) occupation number, 
and (c) 2 component of the total spin vs voltage for voltages of 
amplitude Vb = 0.5 V and relatively short periods T = 1-25, 
2.5, 5, and 10 T -1 . The other parameters are as in Fig. M 
A single hysteresis loop in the periodic regime is shown. The 
arrows indicate the directions in which the loops are traversed. 
The black dashed curves show the voltage dependence in the 
stationary state. 



FIG. 4: (Color online) (a) Current, (b) occupation number, 
and (c) z component of the total spin vs voltage for voltages 
of amplitude Vo = 0.5 V and intermediate periods T = 20, 
40, 80, and 160 The other parameters are as in Fig. |2| 
A single hysteresis loop in the periodic regime is shown. The 
arrows indicate the directions in which the loops are traversed. 
Black dashed curves: Stationary state. 



its initial state. This restriction is relaxed for dynam- 
ical measurements, where the initial state can be tran- 
siently occupied even when the voltage is not sufficiently 
high to populate it in the stationary regime. Thus, ad- 
ditional spectroscopic information on ETE transition en- 
ergies and lifetimes (from the width of the current steps) 
can be obtained from dynamical measurements. 

In the limit of slow driving, T — > 00, the hysteresis 
loops must also close since the observables approach the 
stationary values. The current loops in Fig. |5|a) indeed 
close for slow driving. In particular, the capacitive re- 
sponse (open loop at V — 0) decreases rapidly together 
with the charge at zero voltage. However, the charge and 
the spin in Figs.[5jb) and[5]jc), respectively, have not ap- 



proached the stationary curve even at T — 2560 T 1 . The 
reason for this requires some discussion. 

All quantities show steplike behavior for large T with 
two pairs of steps at voltages of about ±0.20 V and about 
±0.37 V. It is useful to refer to the level scheme in Fig. 
[T] to understand what happens at these voltages. At 
|V| = Vi = 2(E 1>±5/2 - E 0t ± 2 ) = 0.2 V, the excess en- 
ergy eV/2 of electrons appearing in Eqs. (l8])-( 11 ) matches 



the lowest transition energy £^±5/2 — -Eu±2 out of the 
ground states with occupation number n = and mag- 
netic quantum numbers m = ±2 to the states n = 1, 
m = ±5/2. These transitions are denoted by solid ar- 
rows in Fig.[l] From the excited states n = 1, m = ±5/2, 
the molecule can only relax back to the ground states by 
emitting one electron into the leads. No other transitions 
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FIG. 5: (Color online) (a) Current, (b) occupation number, 
and (c) z component of the total spin vs voltage for voltages 
of amplitude Vo — 0.5 V and relatively long periods T = 320, 
640, 1280, and 2560 T~ . The other parameters are as in Fig. 
[2] A single hysteresis loop in the periodic regime is shown. 
The arrows indicate the directions in which the loops are tra- 
versed. Black dashed curves: Stationary state. 



are allowed for sequential tunneling. Thus the molecule 
cannot overcome the anisotropy barrier and the imbal- 
ance AM = (sgn(5 t z t )) between positive and negative 
to cannot be relaxed.^ In fact, this statement is only 
rigorously true at zero temperature. At finite temper- 
atures, there is a thermally activated process in which 
the molecule goes from the ground states to a state with 
n = 1, to = ±3/2 (dashed arrows in Fig.[T]), but the tran- 
sition rate for this process is exponentially suppressed by 
the tail of the Fermi function. 

At the second step at \V\ = V 2 = 2(E l ±z / 2 ~ Eq^± 2 ) — 
0.367763 V, the excess energy eV/2 matches the transi- 
tion energy E 1 ± 3 f 2 — £-o,±2 from the ground states to a 
state with n — 1 and to = ±1/2 (dashed arrows in Fig. 
pj). But then the molecule can overcome the anisotropy 
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FIG. 6: (Color online) Details of the current vs voltage curves 
for periods T = 1.25, . . . , 2560 F -1 from Figs.[3][4] and[H] The 
arrows indicate the direction in which the hysteresis loops are 
traversed. Black dashed curve: Stationary state. 



barrier by a series of sequential-tunneling transitions that 
are either exothermal or endothermal with a transition 
energy lower than eV2/2. Consequently, for voltages close 
to V% , the rate for spin relaxation over the barrier crosses 
over from exponentially small to a sizable fraction of 
rP^ Relaxation over the barrier is still slower than re- 
laxation not crossing the barrier since it involves several 
sequential-tunneling transitions, each of which competes 
with a transition in the opposite direction. 

We can now understand the behavior of the spin for 
slow driving in Fig. [5j Above the second threshold, 
V > V 2 , relaxation over the barrier is not exponentially 
suppressed, and the system indeed approaches the sta- 
tionary behavior for slow driving. As the voltage is low- 
ered into the range Vi < V < V2, the imbalance AM 
between positive and negative to is frozen in. As long 
as the driving is not exponentially slow, the system will 
relax under the constraint of fixed AM. In the range 
Vi < V < V 2 the frozen-in value of AM is close to the 
stationary value so that the system can still relax to a 
state close to the stationary one. When the voltage falls 
below V\, all transitions out of the two ground states 
in Fig. [I] are thermally suppressed, and the system re- 
laxes toward the ground states with the imbalance AM 
approximately conserved. 

For — Vi < V < —Vi, the lowest-energy transitions be- 
come active again but now the frozen value of AM is 
very different from the stationary one at negative volt- 
ages. The system is still with a high probability on the 
left-hand side of the barrier (AM < 0). The most rele- 
vant transition is the one from the state n = 0, to = — 2 to 
the state n — 1, m = —5/2, which requires a spin-down 
electron. However, for negative voltages predominantly 
spin-up electrons are injected. The transition rate for 
this process is thus suppressed by the spin-polarization 
factor p in Eq. ([9]). Consequently, the average occupation 
(n) is suppressed relative to the stationary case, leading 
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to the plateau at (n) ~ 0.336 in Fig. ^b). Finally, for 
V < —V%, spin relaxation over the barrier becomes active 
again, the imbalance AM is unfrozen, and the system can 
relax to the stationary state at slow driving. 

In summary, the device does approach the stationary 
regime in the limit of slow driving, T — > oo, but to reach 
this regime, the period T has to be large compared to the 
exponentially long spin relaxation time. We thus find a 
separation of time scales: The spin relaxation time is 
generically long compared to typical relaxation times for 
processes not crossing the anisotropy barrier. This can 
be compared to the dynamics of a mole cule w ithout local 
spin but involving a vibrational modeP^H j n this case 
small Franck- Condon matrix elements suppressing low- 
energy transitions can lead to a separation of time scales. 



B. Dependence on the voltage amplitude 

Since the imbalance AM between positive and neg- 
ative magnetic quantum numbers is effectively frozen in 
when the voltage drops below the threshold V 2 , we expect 
the behavior of the device to change dramatically when 
the amplitude Vq of the bias V(t) — Vq smut is tuned 
through V 2 . To exhibit this dependence, we have deter- 
mined the periodic behavior for amplitudes Vq through 
the threshold V2, keeping all other parameters fixed. Fig- 
ure [7] shows single hysteresis loops for current, charge, 
and spin in the periodic regime. As expected, the spin 
hysteresis loops are nearly closed for voltage amplitudes 
Vo below the threshold V 2 and open up above the thresh- 
old. Then the system can overcome the barrier for part 
of the driving period so that the spins injected from the 
magnetized leads can reverse the local spin. 

A reasonable measure of the size of the spin hysteresis 
loop — which represents the effectiveness of the device in 
storing information — is the frozen spin at V = with 
V'(t) > 0. We plot the local and the electronic con- 
tributions to the frozen spin as functions of the voltage 
amplitude for zero and finite temperatures in Fig. [HJ To 
obtain the values at zero temperature, we proceed differ- 
ently from the method outlined in Sec. |III| For /cb9 = 0, 
the Fermi function becomes a step function and thus the 
transition rates in Eqs. ([8]) — (11) are piecewise constant 



functions of the instantaneous voltage V(t) and, there- 
fore, of time. This allows us to obtain the stochastic 
matrix II analytically as a product of a finite number of 
matrices describing the time evolution over time intervals 
with constant rates. 

The frozen local spin (S z )q in Fig. |8ja) shows crit- 
ical behavior for vanishing temperature. Taking the 
square we find that the critical exponent is 1/2, (S z )o ~ 
(Vo — V 2 ) x / 2 . The singularity is smeared out at finite 
temperatures. For the frozen electronic spin, critical be- 
havior is not evident. 

The origin of the exponent 1/2 is that the fraction of 
time during which the voltage is large enough to over- 
come the barrier scales with (Vq — V2) 1 / 2 provided that 




-0.2 0.2 

bias voltage V (V) 

FIG. 7: (Color online) (a) Current, (b) occupation number, 
and (c) 2 component of the total spin vs voltage for voltage 
amplitudes V = 0.35V, 0.36V, 0.367763V (= V 2 ), 0.37V, 
0.40 V, and 0.50 V and period T = 40 T -1 . The other pa- 
rameters are as in Figs. [2|j5] A single hysteresis loop in the 
periodic regime is shown in each case; the loop at the thresh- 
old amplitude, Vo = V2, is highlighted as a heavy red (gray) 
curve. The arrows indicate the directions in which the loops 
are traversed. The black dashed curves show the voltage de- 
pendence in the stationary state. 



V(t) is analytic close to its extrema. Taking, for example, 
the first voltage maximum, the end points of this time in- 
terval are obtained by solving V(T/4 ± At/2) = V 2 . If 
V(t) is analytic, we can expand around the maximum, 
V Q + V"(T/A)At 2 /8 = V 2 , the solution of which gives 
At ~ (Vo - V 2 ) 1/2 . We now expand the stochastic ma- 
trix II for small At and use perturbation theory to find 
the probability vector for the periodic state satisfying 
ppcr _ jjP pcr . The leading correction is linear in At 
and thus proportional to (Vq — V 2 ) x / 2 . 

We can therefore draw two conclusions: (i) Quite gen- 
erally, all observables should inherit a (Vq — V2) 1 / 2 cor- 
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FIG. 8: (Color online) (a) Frozen local spin {S z )o and (b) 
frozen electronic spin (s z )o at V(t) = and V'(t) > vs volt- 
age amplitude at thermal energies ksO ~ 0, 2meV, 20meV, 
and 0.2 eV. The other parameters are as in Fig. [7] 

rection from pp° r above the threshold, (ii) The same ar- 
gument also applies whenever additional transitions be- 
come energetically allowed at higher voltage amplitudes. 
Thus there should be corresponding singular terms asso- 
ciated with these transitions. We have checked that this 
is borne out by the numerical results. 




bias- voltage period T (T ) 



spin (S z )o as a function of the voltage period T for three 
temperature values. The voltage amplitude Vq is slightly 
above the threshold V 2 for spin relaxation. At rapid driv- 
ing, T — > 0, the frozen spin goes to zero, as expected since 
the system cannot follow the rapidly changing bias. For 
very large periods, the frozen spin should approach the 
stationary value at zero bias, which is also zero. For the 
(unphysically) high temperature fcs© = 0.2 eV, Fig. [9] 
indeed shows this behavior. On the other hand, at zero 
temperature, the frozen spin never approaches zero for 
large T since there is no spin relaxation for \V(t)\ < V2 
so that the system can never reach the stationary state 
with zero average spin at zero bias. 




bias-voltage amplitude V (V) 



FIG. 10: (Color online) (a) Solid curves: Five periods of the 
z component of the total spin, {S* ot ), vs voltage for differ- 
ent initial pure states with occupation number n = and 
magnetic quantum numbers m = 2, 1, 0, —1, and —2. The 
voltage amplitude is Vo = 0.35 V < V2 and the voltageperiod 
is T = 40 r _1 . The other parameters are as in Figs.[2}j5l Dot- 
ted curve: The same for the initial state with m = but for 
the voltage shifted in time by half a period, (b) Solid curves: 
Five periods of the current for the two cases with m = 2 and 
m = —2. The parameters are the same as in (a). The black 
dashed curves show the true periodic state that the system 
would reach after an exponentially long time. 



FIG. 9: (Color online) Frozen local spin {S z )o vs voltage pe- 
riod at thermal energies fcg9 = 0, 20meV, and 0.2 eV. The 
voltage amplitude is Vo = 0.37 V, slightly above the threshold 
V2 . The other parameters are as in Fig. [7] 

We now turn to the dependence of the frozen spin on 
the driving frequency. Figure [9] shows the frozen local 



C. Quasiperiodic regime 

As we have seen, the anisotropy barrier leads to a sep- 
aration of time scales for relaxation over the barrier and 
relaxation staying on one side of the barrier. It thus 
makes sense to consider the intermediate regime reached 
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after the fast transients have died out, but before the 
slow relaxation has become effective. We call this the 
"quasiperiodic" regime. 

As discussed above, the true periodic state is unique 
since it is determined by the stationary solution of an 
ergodic discrete-time Markov process. Thus the system 
will eventually converge to the same periodic solution re- 
gardless of its initial state. However, the same is not true 
for the quasiperiodic state, which will depend on the ini- 
tial condition and, in principle, also on the protocol with 
which the bias is switched on. This is illustrated in Fig. 



10 'a), which shows the first five periods of the time evo- 
lution of the total spin (S* ot ) starting from different ini- 
tial states. The voltage amplitude is Vq — 0.35 V < V% so 
that the relaxation over the barrier is exponentially slow. 
After a few periods, a quasiperiodic regime is reached, 
which indeed depends on the initial state. After a much 
longer time, during which the system can relax over the 
barrier, all curves would approach the true periodic hys- 
teresis loop (dashed curve). 

We have also verified that the quasiperiodic loop does 
depend on the way the voltage is switched on. The dot- 
ted blue (gray) curve in Fig. [To{a) shows the spin for the 
same initial state as the solid blue (gray) curve; the only 
difference is that the voltage is shifted in time by half a 
period, i.e., V(t) = —Vosinojt. Since the system starts 
in the state with m = 0, which is the one on top of the 
anisotropy barrier in Fig. [T] the sign of the voltage ap- 
plied during the first half period determines the dominant 
spin direction of the injected electrons. This determines 
the probabilities for the system to relax into states with 
positive or negative m. The resulting imbalance AM is 
then frozen in because the voltage amplitude is below the 
threshold Vz. Consequently, the sign of V(t) during the 
first half period determines the spin polarization. 
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FIG. 11: (Color online) Solid red (gray) curve: 50 periods of 
the z component of the total spin vs voltage for the initial pure 
state with occupation number n — and magnetic quantum 
number m = 2. The voltage amplitude is Vo = 0.37 V > V% 
and the voltage period is T = 40F -1 . The other parameters 
are as in Fig. |10| Dashed black curve: Periodic state. 

It is illuminating to compare the behavior when the 



voltage amplitude is above the threshold V%. In this case, 
we do not expect a separation of time scales. Figure [TT] 
indeed shows that the spin approaches the periodic hys- 
teresis loop without reaching any intermediate quasiperi- 
odic regime. The relaxation over the barrier is still slow 
since it involves several sequential-tunneling transitions 
and is only active for a small fraction of the time. 
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FIG. 12: (Color online) z component of the total spin vs (a) 
voltage and (b) time, for a process starting from the pure state 
with n = and m = and consisting of ten periods at the 
amplitude V = 0.35 V < V 2 , five periods at V = 0.5 V > V 2 , 
and five periods at Vo = 0.35 V < V2. The voltage period is 
held fixed to T = 40 T -1 . The other parameters are as in the 
previous figures. 

The spin polarization in the intermediate, quasiperi- 
odic regime can be read out by a transport experiment: 
Figure 10 'b) shows the current hysteresis loops for the 
two cases with initial value m = ±2 (solid curves) and 
also the true periodic behavior (dashed curve). The loops 
are clearly different. To use our model system as a mem- 
ory device, we also need a protocol for writing the spin. 
This is possible by increasing the voltage amplitude over 
the threshold, and making the period T sufficiently long. 
Then the spin approaches a large hysteresis loop on a 
time scale of a few periods. Reducing the voltage ampli- 
tude below the threshold while the spin is large in mag- 
nitude freezes the imbalance AM. For illustration, we 
plot in Fig. [12] the spin for a process starting from the 
pure state with n = and m = and consisting of ten 
periods at a small voltage amplitude Vo < V2, followed by 
five periods at Vq > V2, and eventually by five more pe- 
riods at the smaller amplitude. Figure [12] shows the spin 
as a function of bias and of time. This protocol clearly 
switches the system between two distinct quasiperiodic 
states. If we had reduced the amplitude half of a period 
earlier or later, i.e., when (S^ ot ) was negative, a negative 
spin would have been written. 

The typical time scale of the switching in Fig. [12] is a 
few times T = 40 T -1 - 4 x 10~ 8 s. This should be com- 
pared to the switching time of memristivc systems con- 
taining a MgO layer between ferromagnetic electrodes, 
which is on the order of seconds, probably becau se it 
involves the displacement of oxygen vacancies! 62 * 63 ! The 
switching mechanism proposed in the present paper can 
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be much faster since it involves neither nuclear motion, 
as Refs. [62] and [63] and also the spin transition consid- 
ered by Miyamachi et a/P^ do, nor the motion of domain 
walls, which is the mechanism studied by Miinchenberger 
et oiP 



V. SUMMARY AND CONCLUSIONS 

Memristive (memory resistive) properties of a SMM 
weakly coupled to two ferromagnetic leads have been in- 
vestigated. To that end, we have studied three observ- 
ables: the average current through the molecule, the oc- 
cupation (charge), and the z component of the total spin 
on the molecule. We have obtained the hysteresis loops 
for these quantities vs the applied bias. The device is not 
a purely memristive system, as evidenced by the open 
hysteresis loop for the current. Instead, the transient 
charging of the molecule leads to a partially capacitive 
response. This capacitive response is governed by the 
fast charge relaxation. The device combines memristive 
with spintronics functionality facilitated by a large po- 
larization of the molecular spin. 

For rapid driving, i.e., for driving frequencies on the 
order of the characteristic tunneling rate T, the memory 
dependence is suppressed and the hysteresis loops close. 
The current-voltage characteristics nevertheless show ad- 
ditional spectroscopic features not seen in the stationary 
(dc) current-voltage curve. In the opposite limit of very 
slow driving, memory effects are also suppressed since the 
system eventually approaches the stationary behavior. 
However, the period of the voltage required to reach this 
stationary regime can be exceedingly long in the pres- 
ence of easy-axis anisotropy. The origin of this dynami- 
cal slowing down is that relaxation of the molecular spin 
over the anisotropy barrier becomes suppressed below a 
certain bias voltage so that any imbalance between pos- 
itive and negative spin polarization is essentially frozen 
in. It is the very slow relaxation over the barrier that 
governs the eventual closing of the hysteresis loops. 

At zero temperature, the frozen occupation number 
and spin at zero bias exhibit non-analyticities when the 
voltage amplitude crosses the threshold for relaxation 
over the barrier. The non-analytic contributions scale 
with the voltage distance from the transition point with 
a critical exponent of 1/2. The singular behavior is 
smeared out at finite temperatures. 

We have seen that the easy-axis anisotropy naturally 
leads to a separation of time scales: If the bias voltage 
is below a critical threshold, relaxation over the barrier 
is exponentially suppressed compared to relaxation be- 
tween states on the same side of the barrier, by the tail of 
the Fermi function appearing in the sequential-tunneling 
rates. Now, the question arises as to whether any pro- 
cesses neglected in our approach become relevant in this 
regime. Sequential tunneling is of the order of t^ yh oc T. 
At the order t^ yh cx T 2 , cotunneling occurs: An electron 
can tunnel coherently from one lead to the other. During 



this process, it can flip its spin, which leads to a transition 
of the molecular state without a change of the occupation 
number, but with a change of the spin by one unit. Out 
of the ground states, this process becomes active when 
eV matches the energy difference eV co t = £?o,±i — -Eh, ±2 
between the states with m = ±1 and m = ±2, which 
for the parameters chosen here happens at a voltage be- 
low the Coulomb-blockade threshold V\ , see Fig.[l] (Note 
that the full bias eV enters in this case.) Beyond the volt- 
age Vcat, the system can overcome the anisotropy barrier 
by cotunneling. Thus there is a regime where relaxation 
over the barrier is dominated by cotunneling, which is 
down by a factor of T compared to sequential tunneling, 
but does not involve an exponentially small factor at low 
temperatures. At smaller voltages, | < V co t, cotun- 
neling is also exponentially suppressed and processes of 
even higher order in T become important. Eventually, 
the direct transition from one ground state to the other 
involving a change of the spin by four units occurs at 
order T 8 . This process happens even at zero bias and is 
thus never suppressed by exponential factors. 

Still, below the threshold V 2 for spin relaxation due to 
sequential tunneling, the relaxation rate is at least sup- 
pressed by a factor of F. Thus, for sufficiently weak cou- 
pling between the molecule and the leads, there is still 
a wide separation of time scales. We finally note that 
transverse spin-anisotropy terms, or a transverse mag- 
netic field, can lead to spin tunneling through the barrier 
and thereby open another channel for spin relaxation. 

The separation of time scales opens up an intermedi- 
ate time regime where fast transients have died out, but 
relaxation over the barrier has yet to become effective. 
In this regime, a quasiperiodic dependence of all observ- 
ables on the bias is observed. Unlike the true periodic 
state, the quasiperiodic hysteresis loop depends on the 
initial conditions and the protocol by which the bias is 
switched on. In view of the long lifetime that can be real- 
ized experimentally, this means that this property can be 
used to store information. Indeed, we have demonstrated 
that this spin information can be read out by measuring 
the alternating current, and that it can be rewritten at 
will by judiciously changing the voltage amplitude. Note 
that we have discussed a two-terminal device. None of 
this functionality requires a gate electrode (although the 
latter would add extra flexibility) , which should make the 
practical implementation significantly easier. 

Finally, it is also worth noting that the presence of dif- 
ferent time scales and long relaxation times makes molec- 
ular magnets ideally suited for a host of neuromorphic 
applications,^ ranging from learning circuits^ and as- 
sociative memory^l to the massively parallel solution of 
optimization problems.^ While neuromorphic and mem- 
ristive devices based o n mag netic solid-state structures 
have been sugg 

ested mmm mo lecular magnets poten- 
tially offer higher integration densities, faster switching, 
and lower power consumption. These features make them 
ideal candidates for neuromorphic computing. 
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